
//本法采用了逆Broyden方法

#include <stdio.h>
#include <math.h>
#define E 0.5*1e-6
#define T 1/(6*exp(1)+4)

int main()
{
    double X0[3][1], H0[3][3], r[3][1], F0[3][1], F1[3][1], HF[3][1], y[3][1], Hy[3][1], rH[1][3], rHy[3][1], A;
    double X1[3][1] = {1.0, 1.0, 1.0}, H1[3][3] = {{(exp(1) + 3)*T, (2 * exp(1) - 1)*T, 7 * T}, {(exp(1) + 1)*T, (2 * exp(1) + 1)*T, T}, { -2 * exp(1)*T, 2 * exp(1)*T, 4 * T}};
    int i, j, num = 0;
    do
    {
        for (i = 0; i < 3; i++)             //每次循环，对解矩阵x和矩阵H初始化
        {
            for (j = 0; j < 3; j++)
            {
                H0[i][j] = H1[i][j];
                if (j == 0)
                    X0[i][j] = X1[i][j];
            }
        }

        F0[0][0] = X0[0][0] * X0[1][0] - pow(X0[2][0], 2) - 1;                      //得到F(x)值的矩阵
        F0[1][0] = X0[0][0] * X0[1][0] * X0[2][0] + pow(X0[1][0], 2) - pow(X0[0][0], 2) - 2;
        F0[2][0] = exp(X0[0][0]) + X0[2][0] - exp(X0[1][0]) - 3;

        for (i = 0; i < 3; i++)
        {
            HF[i][0] = 0;
            for (j = 0; j < 3; j++)
            {
                HF[i][0] = HF[i][0] + H0[i][j] * F0[j][0];      //得到矩阵H与F(x)值矩阵的乘积的矩阵
            }
            X1[i][0] = X0[i][0] - HF[i][0];
            r[i][0] = X1[i][0] - X0[i][0];
        }

        F1[0][0] = X1[0][0] * X1[1][0] - pow(X1[2][0], 2) - 1;
        F1[1][0] = X1[0][0] * X1[1][0] * X1[2][0] + pow(X1[1][0], 2) - pow(X1[0][0], 2) - 2; //得到新F(x)值矩阵
        F1[2][0] = exp(X1[0][0]) + X1[2][0] - exp(X1[1][0]) - 3;

        A = 0;
        for (i = 0; i < 3; i++)
        {
            Hy[i][0] = 0;
            rH[0][i] = 0;
            for (j = 0; j < 3; j++)
            {
                y[j][0] = F1[j][0] - F0[j][0];
                Hy[i][0] = Hy[i][0] + H0[i][j] * y[j][0];
                rH[0][i] = rH[0][i] + r[j][0] * H0[j][i];
            }
            rHy[i][0] = r[i][0] - Hy[i][0];
            A = A + rH[0][i] * y[i][0];

        }

        for (i = 0; i < 3; i++)
        {
            for (j = 0; j < 3; j++)
            {
                H1[i][j] = H0[i][j] + rHy[i][0] * rH[0][j] / A;
            }

        }

        num++;

    } while ((fabs(X1[0][0] - X0[0][0]) >= E) || (fabs(X1[1][0]) - X0[1][0]) >= E || (fabs(X1[2][0]) - X0[2][0] >= E));

    printf("The root is:\nx = %.10lf\ny = %.10lf\nz = %.10lf\n\n", X1[0][0], X1[1][0], X1[2][0]);
    printf("The number of iterations is %d.\n", num);
    return 0;
}
